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Abstract 

This paper develops a mathematical model describing the influence that conjugation-mediated 
Horizontal Gene Transfer (HGT) has on the mutation-selection balance in an asexually reproducing 
population of unicellular, prokaryotic organisms. It is assumed that mutation-selection balance is 
reached in the presence of a fixed background concentration of antibiotic, to which the population 
must become resistant in order to survive. We analyze the behavior of the model in the limit of 
low and high antibiotic-induced first-order death rate constants, and find that the highest mean 
fitness is obtained at low rates of bacterial conjugation. As the rate of conjugation crosses a 
threshold, the mean fitness decreases to a minimum, and then rises asymptotically to a limiting 
value as the rate of conjugation becomes infinitely large. However, this limiting value is smaller 
than the mean fitness obtained in the limit of low conjugation rate. This dependence of the mean 
fitness on the conjugation rate is fairly small for the parameter ranges we have considered, and 
disappears as the first-order death rate constant due to the presence of antibiotic approaches zero. 
For large values of the antibiotic death rate constant, we have obtained an analytical solution for 
the behavior of the mean fitness that agrees well with the results of simulations. The results of this 
paper suggest that conjugation-mediated HGT has a slightly deleterious effect on the mean fitness 
of a population at mutation-selection balance. Therefore, we argue that HGT confers a selective 
advantage by allowing for faster adaptation to a new or changing environment. The results of this 
paper are consistent with the observation that HGT can be promoted by environmental stresses 
on a population. 
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I. INTRODUCTION 



Horizontal Gene Transfer (HGT) is considered to be any form of direct transfer of genetic 
material between two organisms, where one organism is not the parent of the other (the 
latter case is known as vertical gene transfer) (Ochman et al. 2000). HGT has become a 
subject of great interest for both molecular and evolutionary biologists, because it is believed 
that HGT plays a large role in re-shaping prokaryotic genomes (Ochman et al. 2000). In 
particular, HGT is believed to be primarily responsible for the rapid spread of antibiotic 
drug resistance in bacterial populations (Walsh 2000). Given that the emergence of antibiotic 
drug resistant strains of bacteria has become a major public health concern (Walsh 2000), 
an understanding of HGT is not only important for advancing our knowledge of biology, but 
it is also of immense practical significance. 

Currently, there are three known mechanisms by which HGT occurs (Ochman et al. 
2000): 

1. Transformation: When an organism (generally a bacterium) collects genetic material 
from its environment. 

2. Transduction: When a virus directly infiltrates a bacterium with genetic material. 

3. Bacterial Conjugation: When a bacterium transfers genetic information via intercel- 
lular contact with another bacterium. 

Bacterial conjugation is believed to be the most important mechanism responsible for 
HGT (Ochman et al. 2000), and so, in this paper, we will focus on developing mathematical 
models describing the role that conjugation- mediated HGT has on the mutation-selection 
balance of bacterial populations. Given the presumed importance that HGT has for the 
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spread of antibiotic drug resistance in bacterial populations, the mathematical models we 
develop will look at the influence of HGT on the mutation-selection balance in the presence 
of an antibiotic. 

The best characterized bacterial conjugation system is the F + /F~ system (Russi et al. 
2008). Here, a bacterium containing what is termed an F-plasmid fuses with a bacterium 
lacking the F-plasmid. The bacterium containing the F-plasmid is termed an F + bacterium 
while the bacterium that does not contain this plasmid is termed an F~ bacterium. When 
the F + bacterium meets an F~ bacterium, it transfers one of the strands of the F-plasmid 
to the F~ bacterium via a pilus. Once a strand of the F-plasmid has been transferred from 
the F + bacterium to the F - bacterium, a copy of the plasmid in both cells is produced by 
daughter strand synthesis using the DNA template strands. The F - bacterium then becomes 
an F + bacterium that transcribes its own pilus and is able to transfer the F + plasmid to 
other bacteria in the population (Russi et al. 2008). This process is illustrated in Figure 1. 

The F + /F~ system is not the most common form of bacterial conjugation. It is what is 
known as a narrow spectrum conjugation mechanism (Tenover 2006), since the F~ plasmid 
may only be transferred between cells that are from similar strains. However, it is known 
that the genes for resistance to various antibiotic drugs have been transferred between dis- 
tinct strains of bacteria, suggesting that a broad spectrum conjugation mechanism is likely 
the important form of HGT leading to the spread of antibiotic drug resistance in bacterial 
populations (Tenover 2006). Nevertheless, because all of the bacterial conjugation mecha- 
nisms follow a pathway that is similar to the F + /F~ pathway, we will use the F + /F~ system 
as the basis for developing our mathematical models of conjugation- mediated HGT. 
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FIG. 1: Illustration of the process of bacterial conjugation. In steps 1 and 2, an F + bacterium 
containing the F-plasmid (blue) binds to an F~ bacterium lacking the plasmid. One of the template 
strands from the F-plasmid then moves into the F~ bacterium, as shown in step 3. In step 4, the 
complementary strands are synthesized to reform the complete F-plasmids in both bacteria. Both 
bacteria are now of the F + type. 

II. MATERIALS AND METHODS 

We assume an asexually reproducing bacterial population, where the genome of each 
bacterium consists of two double-stranded, semiconservatively replicating DNA molecules. 
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The first DNA molecule contains all of the genes necessary for the proper growth and 
reproduction of the bacterium itself. This DNA molecule corresponds to the large, circular 
chromosome that defines the bacterial genome. We assume that there exists a wild-type 
genome characterized by a "master" DNA sequence. It is assumed that a bacterium with 
the master genome has a wild-type fitness, or first-order growth rate constant, given by 1. 
Such a bacterium is termed viable. Furthermore, making what is known as the single-fitness- 
peak approximation (Tannenbaum and Shakhnovich 2005), we assume that any mutation to 
the bacterial genome renders the genome defective, so that the bacterium then has a fitness 
of 0. Bacteria with defective genomes are termed unviable. 

The second DNA molecule is the F-plasmid, which we assume consists of two regions. 
The first region comprises the various genes necessary for bacterial conjugation itself, i.e. 
for allowing the plasmid to move between bacteria. The second region is assumed to encode 
for the various enzymes conferring resistance to a given antibiotic. For this initial study, 
we are interested in the interplay between conjugation- mediated HGT and antibiotic drug 
resistance at mutation-selection balance (we will consider adaptive dynamics later), and so 
this is the simplest model that incorporates these various effects. 

As with the single-fitness-peak approximation made for the bacterial genome, for the F- 
plasmid we assume that there are master sequences for both the conjugation and antibiotic 
drug resistance regions. If the region coding for bacterial conjugation corresponds to a given 
master sequence, then, assuming that the bacterium is also viable, the F-plasmid may move 
into another viable F~ bacterium. Otherwise, we assume that plasmid cannot move into 
another bacterium, in which case the bacterium is treated as an F~ bacterium. 

Similarly, if the region coding for antibiotic drug resistance corresponds to a given master 
sequence, then we assume that the bacterium is resistant to the antibiotic. Otherwise, the 

6 



bacterium is not resistant to the antibiotic, and is assumed to die according to a first-order 
rate constant given by k d . We assume that only viable bacteria interact with the antibiotic, 
since non-viable bacteria do not grow and so may be treated as dead. 

A given genome may be characterized by a three symbol sequence a — ± ± ±, specifying 
the state of the viability, conjugation, and resistance portions of the genome, respectively. 
A "+" is taken to signify that the given genome region is identical to the corresponding 
master sequence, and a "-" is taken to signify that the given genome region differs from the 
corresponding master sequence. 

To develop the evolutionary dynamics equations governing this population, we let n a 
denote the number of organisms in the population with genome a. We wish to develop 
expressions for dn a /dt for the various a. Since we are only interested in the viable population, 
the a of interest are + + +, + H — , H h, H . 

We must now consider the various aspects of the evolutionary dynamics that affect the 
expressions for the dn a /dt. The first aspect of the dynamics that we consider is replication: 
During the semiconservative replication of the bacterial genome, the strands of the DNA 
molecule separate and serve as templates for daughter strand synthesis. Daughter strand 
synthesis is not necessarily error-free, so that there is a probability p, denoted the replication 
fidelity, that a given template strand will produce a daughter genome that is identical to the 
original parent. Because our genome consists of three genome regions, we may define three 
such probabilities, denoted p v , p c , and p r , corresponding to the replication fidelities for the 
viability, conjugation, and resistance portions of the genome. For a replication fidelity p, it 
follows that a template strand derived from a master genome region has a probability p of 
forming a daughter genome region that is identical to the parent, and a probability of 1 — p 
of forming a mutated daughter. If we assume that sequence lengths are long, then making 
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an assumption known as the neglect of backmutations (Tannenbaum and Shakhnovich 2005), 
we assume that a template strand derived from a parent that differs from the master genome 
produces a daughter that differs from the master genome with probability f. The basis for 
this assumption is that for very long genomes, mutations will typically occur in previously 
unmutated regions of the genome, so that mutations will tend to accumulate. 

The second aspect of the dynamics that we consider is conjugation: We assume that 
conjugation occurs between a viable F + -bacterium and a viable F~-bacterium. Thus, con- 
jugation can only occur between a bacterium of type + + ± and a bacterium of type H — ±. 
This process is modeled as a second-order collision reaction with a rate constant 7. The 
conjugation process itself involves the transfer of one of the strands of the plasmid from the 
F + -bacterium to the F~-bacterium, so that the full plasmid needs to be re-synthesized in 
both bacteria via daughter strand synthesis. This introduces the possibility of replication 
errors in either one of the bacteria. 

It should be emphasized that we are assuming for simplicity that all bacteria in the 
population contain exactly one plasmid. This plasmid may contain the correct copies of 
the genes for conjugation, in which case the bacterium is an F + -bacterium, or the plasmid 
may contain defective copies of the genes for conjugation, in which case the bacterium is an 
F~-bacterium. We also assume that, during conjugation, the plasmid transferred from the 
F + -bacterium replaces the plasmid in the F~-bacterium. This is a simplifying assumption 
that will obviously have to be re-examined in future research, where we anticipate developing 
more accurate models that allow for variable plasmid numbers in the bacterial cell. 

Putting everything together, we obtain that the evolutionary dynamics equations are, 



8 



dU jt + = i 2 PvPcPr - 1 + ^(2p c p r - l)(n+_+ + 

= [ZPvPc -1-K D + -^(ZPc - l)(n+-+ + n+__)]n++_ 
+2p c (l -p r )[p„ + — (n + _+ + 

d ~ = [2p«Pr - 1 - y( n +++ + n ++-)) n +-+ + 2(1 - Pc)p r [Pv + y + 

= [2p» - 1 -«r> - + «++-)]«+— + 2(1 -Pc)(l -Pr)[Pv + y( n +-+ +«+--)]«+++ 

7 

+2(1 -Pc)\Pv + ^(ra+_ + + n + __)]n ++ _ + 2p v (l - p r )n+-+ (1) 

where V is defined as the system volume. To put the equations into a form that makes 
the analysis of the mutation-selection balance possible, we define the total population n = 

n +++ + + + n H + + + n + + n , and then define population 

fractions x a via x a = n a /n. We also define a population density p = n/V, and we assume 
that p is constant. Converting from population numbers to population fractions, we obtain, 



dx^ — | — 1_ 



[2p v p c p r - 1 + 1p(2p c p r ~ 1)(X+-+ + - K(t)]x+++ 



dt 

-- [2p v p c -1-k d + lp{2p c - + - re(t)]a; ++ _ 



|_ 



dt 

+2p c (l -p r )\pv +1P(X+-+ +X+— )}X +++ 



dx + _. 



[2p v p r - 1 - lp(x +++ + - K(t)]x+-+ + 2(1 - Pc)p r [Pv + 7P^+- + + )]^+++ 



dt 

[2p v -1-k d - lp(x +++ + x ++ _) - /c(t)]ic + __ 







dt 

+2(1 -Pc)(l -p r )[p„ + 7P(^+-+ + ^+— )]&+++ 

+2(1 -p c )[p„ + 7p(:r + _ + + )]z++- +2p„(l (2) 

where /t(t) = (l/n)(dn/dt) = x +++ +x + ^ + + {\ — ku){x ++ ^+x^ ) is the mean fitness of the 

population. In the subsequent analysis, we will be interested in computing the mean fitness 



at mutation-selection balance, since the mean fitness provides the measure of the effective 
first-order growth constant of the population. Therefore, the mean fitness will allow us to 
understand the selective advantage of HGT in a static environment. 

To determine the values for p v , p c , and p r , we assume that daughter strand synthesis has 
a per-base mismatch probability e, which incorporates all DNA error-correction mechanisms 
such as proofreading and mismatch repair. Because we are assuming complementary double- 
stranded DNA molecules, we assume that all post-replication mismatches are corrected 
via various lesion repair mechanisms (e.g. Nucleotide Excision Repair or NER). However, 
because at this stage there is no discrimination between parent and daughter strands, a 
mismatch is either correctly repaired with probability 1/2, or is fixed as a mutation in the 
genome with probability 1/2. Thus, the net per-base mismatch probability is e/2. If the 
total sequence length is L, then the probability of producing a mutation-free daughter from 
a given parent template strand is (1 — e/2) L . 

If we define \i = Le, so that ji is the average number of mismatches per template strand 
per replication cycle, and if we assume that L — > oo while \i is held constant, then we 
obtain that (1 — e/2) L — > e"^ 2 . For the case of the three-gene model we are considering, we 
let L v , L c , and L r denote the lengths of the genome controlling viability, conjugation, and 
resistance, respectively. Defining L = L v + L c + L r , and a v = L v /L, a c = L c /L, a r = L r /L, 
we then obtain that, 

n — p-otvH/2 
Fv c 

Pc = e -«cM/2 

It should be noted that holding /i constant in the limit of infinite genome length is 
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equivalent to assuming a fixed per genome replication fidelity in the limit of long genomes. 



III. RESULTS AND DISCUSSION 

In this section, we will solve for the mean fitness at mutation-selection balance, denoted 
by R, for two different sets of parameter regimes: We will first consider the case of arbitrary 
kd, but with 7p — > and 7p — > oo. We will then consider the case of arbitrary jp, but 
with kd — > and kd — > oo. Both sets of cases are analytically solvable, and may be used 
to qualitatively understand the behavior of R for arbitrary values of Kd and jp. 

In order to avoid having the derivation of the results interfere with the results themselves, 
for convenience we present the final analytical results for each parameter regime being con- 
sidered, and then provide the derivations in a subsequent subsection. We do not relegate 
the derivations to an appendix, as we believe that they are sufficiently interesting to remain 
part of the main text. 

A. Behavior of R for arbitrary kd 

In the limit where 7p — > 0, the ability for conjugation is lost due to genetic drift (since it 
is never used), and we obtain that, 

k 7P ^ = max{2p„p r - 1, 2p v - 1 - k d } (4) 

We now consider the limit where 7p — > oo. We obtain at steady-state that, 

R _ r lPvPcPr- 1 + 2(1 -Pv)(l-Pc) . 1 

K 7P->oo — maxj — — — ,zp v — 1 — Koj {£>) 

zp c 1 

where x +++ > when R is given by the first expression, and x +++ = when R is given by 
the second expression. 

We can also show that R lp ->oo < ^7p->o- 
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B. Behavior of re for arbitrary 7/) 

Now we consider the behavior of R for arbitrary values of jp, but where Kr> is either very- 
small or very large. Combined with the results of the previous subsection, we may then 
piece together a qualitative sketch of how R depends on Kb and 7p. 

When kd — > 0, there is no selective advantage for maintaining antibiotic drug resistance 
genes in the genome, and so we expect these genes to be lost to genetic drift. Thus, we 
expect, at mutation-selection balance, that x +++ = x + _ + = 0, so we need only consider the 
populations and x H We may also show that R = 2p v — 1. 

Furthermore, the fraction of viable conjugators, x +++ + x ++ -, exhibits a transition as a 
function of 7p. For sufficiently small values of 7p, we have that x +++ + £++_ = 0, while for 
sufficiently large values of 7p, we have that, 

+ x ++ _ = 2 Pv - 1 - 2 ^~-l ^ 
The transition between the two regimes may be shown to occur at, 

< \ 2p v {l- Pc ) 

[1P)trans - (2 Pv -l)(2 Pc -l) [i) 

It may be shown that the disappearance of the conjugators below the critical value of 
7P corresponds to a localization to derealization transition over the portion of the plasmid 
coding for conjugation, so that this transition is a conjugation-mediated HGT analogue of 
the well-known error catastrophe from quasispecies theory (Tannenbaum and Shakhnovich 
2005). 

To understand this behavior, we note that plasmids with defective genes for conjugation 
nevertheless replicate due to the replication of the bacteria in which they reside. Thus, 
for plasmids with functional genes for conjugation to be preserved in the population, their 
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additional growth rate due to conjugation must overcome the loss of functionality due to 
replication mistakes in the genes controlling conjugation. If the conjugation rate is too slow 
and unable to overcome this loss of functionality, then the fraction of conjugators in the 
population drops to zero. 

We now consider the case where kd — > oo. In contrast to the case where 7p — > oo of the 
previous subsection, where we could solve for R for arbitrary values of k d , here we cannot 
readily analytically solve for R for arbitrary values of jp. However, we can obtain analytical 
solutions for R in certain limiting cases of jp, and then interpolate between the two solution 
regimes. As will be seen in the subsection comparing theory and simulation, this approach 
turns out to be fairly accurate. 

In the first limiting case, we assume that 7p remains finite in the limit that kd - ► oo. 
This assures that = 0, since the rate of death due to the presence of antibiotics 

is so fast that no non-resistant genotypes are present in the population. The fact that 7p is 
taken to be finite in the limit that kd — *• oo means that a non-resistant genotype cannot be 
"rescued" via conjugation with a resistant bacterium before death occurs. 

We then obtain that either R = 2p v p r — 1 , or that R is the solution to the following 
equation: 



= 2(1 - p r )R + 2(1 -p v ) (R + l- 2p vPcPr ) 2 

2p c p r - 1 R [1 - 2p r (l - p c )]R - [2p v p c p r - 1 + 2p r (l - p v )(l - p c )] 

In the first case, we have that x +++ = 0, while in the second case we have that x +++ > 0. 
The transition between the two regimes may be shown to occur at, 

, v _ 2p v p r (l - p c )[l - 2p v (l-p r )} 
(2p v p r - l)(2p c p r - 1) 

where x +++ = for 7p < (jp) trans an( i x +++ > f° r IP > (lP)trans- We may show that 
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this expression for (7p)j rans is larger than the corresponding expression for the Kd = case. 

To understand the behavior of R where 7p > (jp) trans, we consider the asymptotic be- 
havior of R in the limit as 7p — > oo. In this case, Eq. (8) reduces to, 

2p v p c p r - 1 + 2p r (l -p v ){l -p c ) f . 

K= ; — yz : (1U) 

l-2p r (l-p c ) V ; 

We may show that this expression is smaller than the expression for R obtained in the 
arbitrary kd, infinite 7p case. 

We now consider the second limiting case in the k d — > oo limit, specifically where 7p is 
itself infinite. Here, however, the ratio between kd and 7p may play an important role in 
the competition between death of non-resistant bacteria, and their "rescue" by conjugation 
with resistant bacteria. Thus, here, we will assume that both 7p, kd — > oo, but we will take 
1P/kd to have some given value in this limit. For large values of this ratio, we expect the 
rescue effect to dominate over bacterial death, and so the value of R should approach the 
value obtained for arbitrary k d in the 7p — > oo limit. For small values of this ratio, we 
expect bacterial death to dominate over conjugation, and so the value of R should decrease 
to a value that will need to be determined. 

We may show that, 

TP = R+_ 2(1 -p v ) [1 - 2p r (l -p c )}R - [2p v p c p r - 1 + 2p r (l - p v )(l - p c )] . . 

K D R [ZPvPcPr - 1 + 2(1 - p v )(l - p c )\ - (2p c - 1)R 

and so obtain that, 

ZPvPcPr - 1 + 2p r (l - p v )(l - p c ) 



l-2p r (l-p c ) 



_ 1p v p c p r - 1 + 2(1 - p v )(l - p c ) 
K 1P/k d -^cg — 2p — 1 

Therefore, for large kd, we expect that R will initially be given by 2p„p r — 1 up to a 
critical value of 7p, after which it begins to decrease according to Eq. (8). Once 7p becomes 
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sufficiently large, we expect that the 7p/kd ratio is such that the functional form for R 
transitions from the finite 7p solution to the infinite 7p, fixed 7p/k_d solution. To estimate 
the transition point between the two solution regimes, we equate the values for 7p as a 
function of R for the two solutions. This allows us to solve for R and thereby allow us to 
solve for 7p. 

We then obtain that the transition point occurs at, 

/ IP y =2p 2 PcPr-l + 2(l-p^)(l-p r ) / p v (l-Pc) , 13 v 

^— trans r 2p v p c p r _ 1 _|_ 2p r (1 - p v ) (1 - p c ) \j 1 - 2p r (1 - p c ) 

Note that, as kd — > oo, we have that (jp) trans ~^ 00 an d {iPl k d) trans —> 0, so the 
assumptions that allowed us to make the calculation above are valid. 

C. Comparison of Theory and Simulation 

Figure 2 shows plots of R versus p for both the 7p — > 0, 7p — > 00 limits. Plots were 
obtained using both the analytical formulas obtained in this paper, as well as via stochastic 
simulations of replicating organisms. Note the good agreement between theory and simula- 
tion. 

Figure 3 illustrates the regimes, as a function of p and 7p, where there exist a positive 
fraction of conjugators at steady-state, and where the fraction of conjugators is zero. This 
is computed for the kd = limit. Note that, as p increases, 7p must be pushed to higher 
values so that there is a positive fraction of conjugators at steady-state. As explained before, 
this increase in 7p is necessary to overcome the mutation-induced loss of functionality as p 
increases. 

Figure 4 shows three plots of R versus 7p for k d = 10. One of the plots was obtained by 
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FIG. 2: Plots of R versus /i for both the 7/5 — > 0, 7/9 — > oo limits. The parameter values we 
took are a v = 0.6, a c = a r = 0.2, and ke> = 10. We show both analytical results and results from 
stochastic simulations. The analytical results are plotted using thin solid lines, where the top curve 
corresponds to the 7/3 = result, while the bottom curve corresponds to the 7/5 = oo result. The 
dotted line corresponds to the stochastic simulation for 7/) = 0, and the dashed line corresponds 
to the stochastic simulation for 7/9 = 00. Parameter values for the stochastic simulations were 
L v = 30, L c = L r = 10, and a population size of 1, 000. 

numerically solving for the mutation-selection balance using fixed-point iteration. The other 
two plots correspond to the infinite kd, finite 7/9, and infinite Kd, fixed 'Jp/kd expressions 
for R given in the preceding subsections. Note that already for k,£> = 10 the approximate 
analytical solutions capture the dependence of R on 7/) fairly accurately. 
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FIG. 3: Regimes of existence and non-existence of conjugators as a function of /x and jp, where 
kd = 0. The boundary between the two regimes was computed analytically. 

D. Derivation Details of the Analytical Results 

1. Derivation of R for arbitrary kd, and 7/3 — > 

Due to the nature of exponential growth, for the population fractions to converge to a 
stable steady-state we must have that, R > 2p v p c p r — 1, 2p v p c — 1 — kd, 2p v p r — 1, 2p v — 1 — kd- 
Because 2p v p c p r — 1 < 2p v p r — 1, and 2p v p c — 1 — Kd < 2p v — 1 — kd, it follows that 
k — %PvPr ~ l ; 2p„ — 1 — Kd. However, if we then look at the steady-state version of Eq. 
(2), obtained by setting the time derivatives to 0, we then obtain that x +++ = a; ++ _ = 0. If 
> 0, then the third equation gives us that R = 2p v p r — 1, otherwise the fourth equation 
gives us R = 2p v — 1 — kd- 
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FIG. 4: Plots of k versus jp for kd = 10, fi = 0.4, a v = 0.6, a c = ct r = 0.2. The plot marked with 
the solid line was obtained by numerically solving for R using fixed-point iteration. The dashed 
line was obtained by using the infinite kd, finite ^p expression for R, while the dotted line was 
obtained by using the infinite kd, fixed jp/kd expression for R. 



So, we have shown that R > 2p v p r — 1, 2p v — 1 — Kd, and yet R = 2p v p r — 1 or 2p v — 1 — kd- 
These two requirements imply that R = msix{2p v p r — l,2p v — 1 — k^}. Note that we have 
also shown that = 0, so that our claim that conjugation is lost due to genetic 

drift has also been proven. 
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2. Derivation of k for arbitrary kd, and 7/) — > oo 

In the limit where 7p — > oo, we have that X _| [_ — X_| = 0. However, 7px + _ + and 

7px H may converge to positive values. So, we define = jpx + _ + and = 7px^ 

Because x + - + = x H = 0, we also have that dx + - + /dt = dx^ /dt = 0, and so from 

Eq. (2) we have that, 

= — z + - + {x +++ + £++_) 

+2(1 -p c )\pv + z+-+ + z + —]p r x++ + 

= + x ++ _) 

+2(1 - p c ) [p v + + + [(1 - Pr)x +++ + 

(14) 

Summing these two equations and solving for 2+_ + + z A gives, 

2(1 -p c )p v 

Z+ -+ + Z+ - = 2p e -l (15) 

Substituting into the expressions for dx +++ /dt and dx ++ _/dt from Eq. (2) we obtain, 

after some manipulation, 

cfe +++ ,2p v p c p r - 1 + 2(1 -p c ) _ M1 

= [ 2^ 

<fo++- ro n 2p„p c (l-p r ) 
rft = l 2 Pv -1-kd- re(t)]a;++_ + — 2 p -\ — X+++ 

(16) 

Following a similar argument to the 7p — > case, we obtain the expression for 
given above. 

To prove o, we need only show that, 

ZPvPcPr - 1 + 2(1 - p v )(l - p c ) 



2p c -l 
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< 2p vPr - 1 (17) 



After some manipulation, it may be shown that this inequality is equivalent to, p r < 1, 
which clearly holds, thereby proving the claim. 

3. Derivation of R for kd — > 0, and arbitrary 7/) 

We can add the first two equations from Eq. (2), and also the third and fourth equations, 
to obtain the pair of equations, 



d ( x + ++ + x++ -^ = [ 2pvPc - 1 + lp (2p c - l)(x + _ + + - R(t)](x +++ + x ++ _) 

^±z±±^±=j = [ 2 p v - 1 - 7p(x +++ + - «(*)](&+-+ + 

+2(1 -p c )[p„ + 7p(x + _+ + + (18) 

Summing these two equations then gives, 



^±±±±f±±=±^±=±±^) = [2j , w _ ! _ m{x+++ + + + (19) 

from which it follows that R = 2p v — 1 at steady-state. 

Substituting this value for R into the steady-state version of Eq. (18), we obtain, 

= [(2p c - l) 7 p(a; + _ + + - 2p v (l - Pc )](x +++ + x ++ _) (20) 

which gives either that x +++ + x ++ _ = or x + ^ + + x A = 2p^(l — p c ) /[ipi^Pc — !)]■ If the 

second case holds, then since 2p v — 1 = R = x +++ + + £+_+ + x H , we obtain that, 

x ++++ x ++ _ = 2p„-l-Ml^ (21) 

Now, for large values of jp, we expect that the population will consist of a non-zero fraction 
of conjugators, so that x +++ +x ++ ^ > 0. However, because cannot be negative, 
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we must have that, 



\ / \ _ 2p v (l-p c ) 
IP > (ip) trans = (2 ^_ 1)(2pc _ 1) (22) 



in order for x +++ + x ++ _ > 0. Therefore, by continuity, we expect that x +++ + x ++ _ = 
for 7p < (7p)< ra „ s , and x +++ + x ++ _ = 2p v - 1 - ^pl-i) > for lf> > Mtrow 



^. Derivation of k for kd — ► oo, and /mite 7p 

In this limiting case, although = 0, it is possible that y++- = k d x ++ - and 

y H = K £) a; H have non-zero, finite values in the limit as k d — > oo, and so we need to 

consider the effect of these quantities in our analysis. We then have that the steady-state 
version of Eq. (2) reads, 



= [2p v p c p r - 1 + 'Jp(2p c p r - - K\X +++ 

= [2p v p r - 1 - ipx +++ - k}x + _ + + 2(1 - p c )p r [p v + -fpx + _ + ]x +++ 

V++- = 2 Pc(l - Pr)[Pv + -1pX + _ + ]x +++ 

V+— = 2(1 -Pc)(l -Pr)[p« + 7P^+-+]^+++ + 2p„(l -p r )x + - + (23) 

If = at steady-state, then R = 2p v p r — 1. So, let us consider the case where 

x +++ > 0. Summing the first two equations from Eq. (23) gives, 

2(1 - p r )-ipx +++ x + _+ = [2p v p r - 1 - «] (x +++ + x + _+) (24) 

Summing the last two equations from Eq. (23) then gives, 

y++- + y+— = [2Pv - i - k](x +++ + x + _ + ) (25) 
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Now, in the limiting case being considered here, we have that k = x +++ + — y++- — 
y+— = [k + 2(1 - p v )\(x +++ + and so, 

X+++ + X+ - += -. + 2(l- Pv ) (26) 

Since x +++ > 0, the first equation from Eq. (23) gives, 

R + 1 - 2p v p c p r 

x +-+ = 7^ 7T~ 27 



and so, 



x = K K + 1 ~ 2p v p c p r 

X+++ ~ R + 2(1 - p v ) 1P (2p cPr -l) 1 ' 



Substituting into Eq. (24) gives the following non-linear equation that k must satisfy: 

9n _ v R + 1 ~ 2p v p c p r R R + l - 2p v p c p r _ R _ _ 

[ Pr) 2p cPr -l [ R + 2(l-p v ) 1P (2p cPr -l) l ~ R + 2(l-p v ) [ZPvPr Kl 

(29) 

which, after some manipulation, may be shown to be equivalent to Eq. (8). 

To determine the critical value for the transition between the x +++ = and x +++ > 
regimes, we note that if x +++ is continuous at this transition, then we must have that 
x +++ = using the expression in Eq. (28), which gives that R = 2p v p r — 1 from Eq. (29), 
so that R is also continuous at this transition. Solving for the critical value of 7p then gives, 

2p v p r {\ -p c )[l - 2p v {\ -p r )\ 

{IP) trans = 77. TT~77\ 7~\ ^ 

(2p v p r - l)[2p c p r - 1) 

So, for 7p < (jp)trans, w e have that x +++ = and R = 2p v p r - 1, while for 7p > (7p) tmns 

we have that x +++ > and R is given by the solution to Eq. (8) or, equivalently, Eq. (29). 

To show that this value for (jp) tTans is larger than the corresponding value obtained for 

Kb = 0, we need to show that, 

2p v p r (l - p c )[l - 2p v (l - p r )} 2p v (l-p c ) 

(2p vPr - l)(2p cPr - 1) "~ > (2p v - l)(2p c - 1) 1 } 
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After some manipulation, this inequality may be shown to be equivalent to, 



Ap v p r (2p c - 1)(1 - p v ) + 2p v p r - 1 > 



(32) 



which clearly holds, and so the inequality is established. 

Finally, to show that the value of R as 7p — > oo is smaller than the value of R obtained 
in the arbitrary kd, 7p — > oo limit, we need to show that, 



which establishes the inequality. 

5. Derivation of R for kd — ► oo, and fixed value of ^p/np 

Because 7p is infinite, we expect that x + _ + = x A = 0, although z + _ + = ^px + _ + 

and z A = 7px H may converge to positive, though finite, values. Also, because the 

+ H — genomes, as conjugators, cannot be "rescued" by conjugators themselves, we expect 
that = in the limit as ke> — > oo, though again it is possible that = K£>x ++ - 

converges to a positive value. We only expect x +++ > 0, since the + + + genomes are both 
conjugators and resistant to the antibiotic, and so are not destroyed by conjugation or by 
antibiotic-induced death. 

The steady-state equations then become, 



tyvPcPr - 1 + 2p r (l - p v )(l - p c ) 
1 - 2p r (l -p c ) 

tyvPcPr - 1 + 2(1 - p v )(l - p c ) 

2p c -l 



(33) 



After some manipulation, this condition may be shown to be equivalent to, 



p v (2p c p r - 1)(1 -p c )(l ~Pr) > 



(34) 
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R = 2p v p c p r - 1 + (2p c p r - l)(z+-+ + 2+—) 
V++- = 2p c (l - p r )[Pt, + Z+-+ + Z+—]X +++ 
Z+-+ = 2(1 - p c )p r [p v + + 

— = [2(1 - p c ){\ - p r )(p v + z + _ + + z+__) - (35) 

IP 

From the first equation we have that + z^ = (« + 1 — 2 PvPcPr ) / (2 PcPr — 1). We 

therefore have that, 



2p c (l -p r ) 

2p c p r - 1 



= 2(l- Pc ) Pr ( _ + i 
2p c p r - 1 

[1 - 2p r (l - p c )]R - [2p v p c p r . - 1 + 2p r (l - p v )(l - p c 



z +-+ = ~7^rz —\ K + 1 ~Pv 



2p c p r - 1 

k d 2p v p c p r - 1 + 2(1 - p v )(l - p c ) - (2p c - l)R 

Z+— = z : x +++ (6b) 

IP 2p c p r - 1 

and we also have in this limit that R = x +++ — y++- — kd/{ip)z^ Substituting in the 

expressions for y ++ - and kd/('Jp)z + , we obtain, 

" +++ = R + 2(1- Pv ) (37) 
Substituting this expression into the last equality of Eq. (36), and using the expression for 
z H , gives us Eq. (11). 



6. Derivation of the transition point between the two functional forms for k for kd — ► oo 

Equating the finite 7p with the infinite 7p expressions for R, we obtain that the transition 
point occurs where, 
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[1 - 2p r (l -Pc)]k - [2p v p c p r - 1 + 2p r (l-p v )(l -p c )] 



K+l- 2p v p c p r 



X 



2(1 -Pr) 

2p c p r - 1 



([2p vPcPr - 1 + 2(1 -p v )(l -p c )\ - (2p c - 1)k) 



(38) 



Since /tp — > oo, we then obtain that the transition point occurs where the left-hand side 
is zero, so that R — [2p v p c p r — l + 2p r (l—p v )(l—p c )]/[l — 2p r (l—p c )]. To estimate the value 
of 7p where this transition occurs in the limit of large k d , we substitute the expression for 
[1 - 2p r (l - p c )\R - [2p v p c p r - 1 + 2p r (l - p v ){l - p c )\ given in Eq. (38) into Eq. (8), and 
then substitute the value of R that we obtained for the transition. After some manipulation, 
we obtain the expression given by Eq. (13). 

IV. CONCLUSIONS 

We have developed a mathematical model describing the role that conjugation- mediated 
Horizontal Gene Transfer (HGT) has on the mutation-selection balance of a unicellular, 
asexually reproducing, prokaryotic population. Because HGT is believed to play a major 
role in the spread of antibiotic drug resistance in bacteria, we considered the effect of an 
antibiotic on the mutation-selection balance of the population. Interestingly, we found that, 
in a static environment at mutation-selection balance, conjugation actually reduces the mean 
fitness of the population. However, by studying the dependence of the mean fitness on 7p 
for large values of k d , the antibiotic-induced first-order death rate constant, we find that 
the behavior is somewhat more complicated: For small values of 7p, the mean fitness is 
constant, and the fraction of viable conjugators in the population is 0. At a critical value of 
7P, the fraction of viable conjugators begins to increase, and the mean fitness decreases to 
its minimum value. After reaching its minimum, the mean fitness increases asymptotically 
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to the 7p — > oo limit, which is nevertheless smaller than the small 7p value for the mean 
fitness. We developed approximate analytical solutions for the functional dependence of 
the mean fitness on 7p in the limit of large kd, and found that these solutions agree well 
with simulation. It is important to note that the fitness variations as a function of 7p were 
fairly small for the parameter values studied. Nevertheless, we believe that this is non-trivial 
behavior that is important to characterize. 

Although the results of our paper are based on a highly simplified model, they nevertheless 
suggest that HGT does not provide a selective advantage in a static environment. This is 
likely due to the fact that, due to mutation, HGT can destroy antibiotic drug resistance 
in a previously resistant cell. While HGT can also confer resistance to a non-resistant cell, 
natural selection alone is sufficient to maximize the population mean fitness in a static 
environment. HGT simply has the net effect of destroying favorable genes, thereby lowering 
the mean fitness. This result may be viewed as an example of the "If it is not broken, do 
not fix it" principle. 

Thus, based on the results of this paper, we argue that HGT likely only has a selective 
advantage in dynamic environments, where it would act to speed up rates of adaptation. 
While this result needs to be checked in future research, it is nevertheless consistent with the 
observation that bacteria can regulate their rates of HGT. For example, it is known that, 
in response to stress, bacteria can activate the SOS response (Beaber et al. 2004), which 
has the effect of increasing rates of HGT. This is consistent with our results suggesting that 
HGT should be kept at a minimal level in static environments, and increased in dynamic 
environments. It is also worth mentioning that while conjugation-mediated HGT has not 
been specifically modeled before in this manner (at least to our knowledge), other HGT-like 
models have been studied (Park and Deem 2007; Cohen et al. 2005), and have found that 
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HGT does indeed allow for faster adaptation in dynamic environments (Cohen et al. 2005). 
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